/* mex 文件对应的 C 代码针对如下运算:
clear;clc;

% N1=30; N2=60; N3=20; M=2000;
N1=30; N2=60; N3=2000; M=2000;

R3 = 30*rand(N3, M);
R2 = 20*rand(N2, M);
R1 = 10*rand(N1, M);
RM = rand(M, 1);

tic
RN=zeros(N1*N2*N3,1);
for k=1:M
temp=kron(R3(:,k),R2(:,k));
RN=RN+kron(temp,R1(:,k))*RM(k);
end
toc

% ------------------------------ 测试
R1=R1.'; R2=R2.'; R3=R3.';

tic
RN_mex=RN_RM123(RM,R1,R2,R3);
toc

% ------------------------------ 比较
isequal(RN, RN_mex) % 有可能因为计算误差而不相等
sum(abs(RN - RN_mex))
*/

#include "mex.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    size_t M, N1, N2, N3;
	double *RN,
		*RM, *RM0, *RM00,
		*R1, *R10, *R1000,
		*R2, *R20, *R200, *R2000,
		*R3, *R30, *R300, *R3000;

	////////////////////////////////// 输入参数检查及预处理 ////////////////////////////////////
	if (nrhs < 4) mexErrMsgTxt("至少 4 个参数!");

	M = mxGetNumberOfElements(prhs[0]); // RM 的元素个数

	if (M != mxGetM(prhs[1]) || // 行数
		M != mxGetM(prhs[2]) ||
		M != mxGetM(prhs[3]))
		mexErrMsgTxt("参数的行数不匹配!");

	N1 = mxGetN(prhs[1]);// 列数
	N2 = mxGetN(prhs[2]);
	N3 = mxGetN(prhs[3]);

	RM0 = mxGetPr(prhs[0]);// 输入矩阵
	R10 = mxGetPr(prhs[1]);
	R20 = mxGetPr(prhs[2]);
	R30 = mxGetPr(prhs[3]);

	//////////////////////////////////////// 输出参数 //////////////////////////////////////////
	plhs[0] = mxCreateDoubleMatrix(N1*N2*N3, 1, mxREAL);

	RN = mxGetPr(plhs[0]);// 输出矩阵

	/////////////////////////////////////// 计算过程 ///////////////////////////////////////////
	R1000 = R10 + M*N1;
	R2000 = R20 + M*N2;
	R3000 = R30 + M*N3;

	RM00 = RM0 + M;

	for (R300 = R30; R300 < R3000; R300 += M)//因为数据是按照列读取的，所以每一行相隔为一个的数据，它们的索引相隔M
		for (R200 = R20; R200 < R2000; R200 += M)//同理，每次读取R2矩阵的特定行的元素进行kron()运算.
			for (R1 = R10; R1 < R1000; RN++)//R1是按照列读取的
				for (R2 = R200, R3 = R300, RM = RM0; RM < RM00;)
					// 注意: 变量相乘顺序不同, 会导致不同的计算误差
					*RN += (*R3++)*(*R2++)*(*R1++)*(*RM++);
}